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Self-organization of charged particles on a 2D lattice, subject to an anisotropic Jahn- Teller-type 
interaction and 3D Coulomb repulsion is investigated. In the mean-field approximation without 
Coulomb interaction, the system displays a phase transition of first order. In the presence of 
the Coulomb repulsion the global phase separation becomes unfavorable and the system shows a 
mesoscopic phase separation, where the size of the charged regions is determined by the competition 
between the ordering energy and the Coulomb energy. 

The phase diagram of the system as a function of particle density and temperature is obtained 
by systematic Monte Carlo simulations. With decreasing temperature a crossover from a disordered 
state to a state composed from mesoscopic charged clusters is observed. In the phase separated 
state charged clusters with even number of particles are more stable than those with odd number of 
particles in a large range of particle densities. With increasing particle density at low temperatures 
a series of crossovers between states with different cluster sizes is observed. Above half filling in 
addition to the low temperature clustering another higher temperature scale, which corresponds to 
orbital ordering of particles, appears. 

We suggest that the diverse functional behaviour - including superconductivity - observed in 
transition metal oxides can be thought to arise from the self-organization of this type. 

PACS numbers: 



I. INTRODUCTION 

The presence of nanoscale inhomogeneities is ubiqui- 
tous in the cuprate superconductors [ a , [| L [a, [a H , the 
magnetoresistive manganites d ©, [Tajnjjl|jl3 and 
other doped transition metal oxidesjliT 15 fig] . Fur- 
thermore, there is emerging consensus that doped charge 
carriers in the oxides may phase segregate to form nano- 
scale textures. These are believed to be of importance for 
achieving their functional properties such is superconduc- 
tivity in the cupratesQ and giant magnetoresistance in 
the manganites [17]. 

For the cuprates the idea of charge segregation ap- 
peared soon after the discovery of superconductivity 
[3 [n| [2(3]. In a doped semiconductor the phase sep- 
aration may have two different origins. The first is the 
chemical origin and is associated with the segregation of 
dopant atoms. This type of phase segregation is usu- 
ally temperature independent and weakly dependent on 
external perturbation. Exceptions may appear due to a 
large mobility of dopant atoms at relatively high temper- 
ature. 

If the mobility of impurity atoms is small one might 
expect a pure electronic mechanism of phase separation. 
In this case the electronic system is in thermodynamic 
equilibrium and competing phases are close in energy. 
This is typical for the systems exhibiting a first order 
phase transition. Electronic phase separation is very of- 
ten observed in the magnetic semiconductors like EuSe 
or EuTe [2 ill 2 21 . 23] . Therefore the idea of the charge seg- 
regation in the cuprate superconductors and in the man- 
ganites is very often associated with magnetic degrees of 



freedom [24|, |25j, |26[ , where the phase separation is dis- 
cussed within t — J model. In Refs. [23, [28|, [29|, [3(| phase 
separation was studied within Hubbard model. The re- 
sults are still controversial. In some cases the t — J model 
displays clear static (24| or dynamic (25[ phase separation. 
The situation is quite different for the Hubbard model. 
For example the results of numerical simulations [30j sug- 
gest that the phase separation is absent at any set of pa- 
rameters and for any size of the lattice. Nevertheless, all 
these models do not consider long-range Coulomb repul- 
sion which has very str ong effect on the phase separation. 

[II [Mill EH 11 If 

The long-range Coulomb repulsion together with the 
surface energy determine the topology of the two phase 
state. The charged carriers have the tendency towards 
spatial segregation which is caused by the fact that the 
free energy density of the phase with finite density of 
carriers is lower then the free energy density of the un- 
doped system. On the other hand, the charge segregation 
leads to the charging effect because the dopant atoms 
are distributed uniformly in the system. Therefore, a 
strong electric field appears which has tendency to mix 
the charged phases. In the low doping limit there is a low 
concentration of charged droplets and they do not over- 
lap. The system behaves as an insulator. When the con- 
centration increases the percolative transition to a new 
phase is expected (17L 1371 138| . 

More recently it was suggested that an interplay 
of a short range lattice attraction and the long-range 
Coulomb repulsion could lead to the formation of short 
metallic or insulating strings of polarons[33. l40j . This 
was mainly motivated by observation of giant isotope ef- 
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feet in manganites and cupratesfUl |42|. In ref. [43| we 
suggested that an anisotropic mesoscopic Jahn- Teller in- 
teraction between electrons and k 7^ optical phonons 
might lead to the formation of carrier pairs and stripes. 
A slightly different approach, based on elasticity was con- 
sidered more recently for the case of the manganites by 
Kugel and Khomskii [3] using the methods of Eremin et 
al. [45| , and by Shenoy et al. [46[ . 

The fundamental question which we try and answer 
here is how charged particles order in the presence 
of anisotropic Jahn- Teller type interaction, particularly 
when their density becomes large. We consider charged 
particles on a 2D square lattice subject to only the 
long-range Coulomb interaction and an anisotropic Jahn- 
Teller (JT) deformation. In the preliminary report we 
have considered a narrow doping range, but have found 
a clear evidence of phase segregation and preferential for- 
mation of pairs. 33] Here we extend this study over the 
full doping range. 

In the mean field (MF) approximation without 
Coulomb repulsion, the system displays a first order 
phase transition to an ordered state below some criti- 
cal temperature. In the presence of Coulomb repulsion 
global phase separation becomes unfavorable and the sys- 
tem shows a mesoscopic phase separation, where the size 
of charged regions is determined by the competition be- 
tween the ordering energy and the Coulomb energy. Us- 
ing Monte-Carlo (MC) simulations we show that the sys- 
tem can form many different mesoscopic textures, such 
as clusters and stripes, depending only on the magnitude 
of the Coulomb repulsion compared to the anisotropic 
lattice attraction and the density of charged particles. 
Surprisingly, in agreement with previous report a feature 
arising from the anisotropy introduced by the Jahn- Teller 
interaction is that in a wide part of the phase diagram 
objects with even number of particles are found to be 
more stable than with odd number particles, which could 
be significant for superconductivity when tunnelling is 
included ■ 

II. FORMULATION 

The model proposed in the ref. [43| involves all interac- 
tions allowed by the symmetry. We consider a simplified 
version of the model, where only the interaction leading 
to the deformation of the B\ g symmetry is taken into 
account. The interaction with Bi g mode leads to simi- 
lar effects and therefore for our purposes we can restrict 
ourselves by consideration B\ g mode only. As a result 
the interacting part of the Hamiltonian has the form: 

Hjt = ff5>s,i{(r* - r 2 y )f Q (r)}(bt +r + h+r), (1) 

r,l 

here the Pauli matrix 03.1 describes two components of 
the electronic doublet, and fo(r) is a symmetric func- 
tion describing the range of the interaction. We omit the 
spin index in the sum, since we ignore spin structure at 



present. The resulting model could be easily reduced to 
a lattice gas model. This is performed using the Lang- 
Firsov transformation or equivalently the adiabatic ap- 
proximation for the phonon field. Let us introduce the 
classical variable $i = + 6j)/v2 and minimize the 
energy as a function of $i in presence of the harmonic 
term $ 2 /2. We obtain the deformation, which cor- 
responds to the minimum energy, 

$. (0) =-V2 5 /u,5>3, £+r /(r), (2) 

r 

where /(r) = (r 2 — r^)/o(r). Substituting $[ ^ to the 
Hamiltonian (1) and taking into account that the car- 
riers are charged we arrive to the lattice gas model. To 
formulate the model we use a pseudospin operator S with 
S = 1 to describe the occupancies of the two electronic 
levels ni and 712. Here S z — 1 corresponds to the state 
with ni = 1 , 712 = 0, Sf = — 1 to ni = 0, TI2 = 1 and 
Sf = to ni = 112 — 0. Simultaneous occupancy of 
both levels is excluded due to the high onsite Coulomb 
repulsion (CR) energy. The Hamiltonian in terms of the 
pseudospin operator is given by 

H JT-c = J2(Vi(i-mS! + V c (i-3)Q l Q i ), (3) 

where Q\ = (S?) 2 . V c (m) = e 2 /eoam is the Coulomb 
potential, e is the charge of electron, eo is the static di- 
electric constant and a is the effective unit cell period. 
The anisotropic short range attraction potential is given 
by 

^(m)= 5 2 /^EiW( m + i )- (4) 

i 

The attraction in this model is generated by the in- 
teraction of electrons with optical phonons. The radius 
of the attraction force is determined by the radius of the 
electron-phonon interaction and the dispersion of the op- 
tical phonons [ioj. 

A similar model can be formulated in the limit of the 
continuous media. In this case the deformation is charac- 
terized by components of the strain tensor. For the two 
dimensional case we can define 3 components of the strain 
tensor: ei = u xx + u yy transforming as the A\ g represen- 
tation of the D^h group, e = u xx — u yy transforming as 
the B\g representation and e% = u xy transforming as the 
Bi g representation. These components of the tensor are 
coupled linearly with the two-fold degenerate electronic 
state which transforms as the E g or E u representation 
of the point group. Similarly to the case of interaction 
with optical phonons we will keep the interaction with 
deformation of the B\ g symmetry, namely e only. The 
Hamiltonian without the Coulomb term has the form: 

H=~gJ2 S "« + \ ( A <i + + Mi) ( 5 ) 
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here Aj are corresponding components of the elastic mod- 
ulus tensor, and g is coupling constant of the charge car- 
riers with the strain tensor. The components of the strain 
tensor are not independent [46[ and obey the compatibil- 
ity condition: 

V 2 ei(r) - Ad 2 e 3 (r)/dxdy = (d 2 /dx 2 - d 2 /dy 2 )e(r). 

The compatibility condition leads to the long range 
anisotropic interaction between polarons. To derive the 
Hamiltonian we minimize Eq.(5) with respect to e\ and 
e 3 taking into account compatibility condition. The re- 
sulting Hamiltonian in the reciprocal space has the form: 



H 



(A 2 + A 1 U(k)) 



(6) 



The wavevector dependence of the potential is given by 



8(A 1 /A 3 )k 2 x kl 



By minimizing the energy with respect to tk 
and including the long-range CR we again ob- 
tain Eq.(3). The anisotropic interaction potential 

Vi(m) = - J2 k exp (ik ■ m) - 



A^i'^-'V ^Affik)) is determined 
in this case by the interaction with the classical defor- 
mation and is long-range as well. It decays as 1/r 2 at 
large distances. Since attraction forces decay faster then 
the Coulomb repulsion at large distances the attraction 
can overcome the Coulomb repulsion at short distances 
leading to the mesoscopic phase separation. 

Irrespective of whether the resulting interaction be- 
tween polarons is generated by acoustic or optical 
phonons the main physical picture remains the same. 
In both cases there is an anisotropic attraction between 
polarons over short distances. This interaction can be 
either ferromagnetic or antiferromagnetic in terms of 
the pseudospin operators depending on the spatial direc- 
tion. Without loosing generality we assume that V(m) is 
nonzero only for the nearest neighbors and can be either 
ferromagnetic or antiferromagnetic. 



III. MEAN FIELD 

Our main goal is to study this lattice gas model (3) at 
a constant average density, 



4i> 



(S) 



where N is the total number of sites. However, to clarify 
the physical picture we first perform calculations in ab- 
sence of long-range CR at a fixed chemical potential first 
by adding the term — /U^i Qi to the Hamiltonian Q. 

Similar models were studied many years ago on the 
basis of the molecular-field approximation in the Bragg- 
Williams formalism [I?], |48| . The mean-field equations 



for the two variables n and M 
form [H: 



7f J2i &i nave tne 



M = 



2sinh (2zViM/k B T) 



cxp {-n/k B T) + 2 cosh (2zV l M/k B T) ' 



2 cosh (2zViM/k B T) 



exp (-n/k B T) + 2 cosh (2zViM/k B T) 



(9) 



(10) 



here z = 4 is the number of the nearest neighbours for 
the square lattice in 2D and k B is the Boltzman constant. 
For positive fj, > equation ((9]) has 2 solutions bellow T c . 
The solution with M = is unstable while the solution 
with a finite M corresponds to the global minimum with 
n — > 1 for T — > 0. When —2zVi < \i < the equation 
has 3 solutions below T c \ < T c . The free energy has two 
minima and one maximum. The phase transition at T c \ 
is of first order. The trivial solution M = corresponds 
to the case when n — > as T — ► 0. For u < — 2zVi there 
is only the trivial solution of the equation M = 0. 

When the number of particles is fixed (Eq.©) 
the system is unstable with respect to global phase 
separation below T cr it(n). The line of the phase 
transition is determined by the condition: F(M = 
Q,fi crit (T),T) = F(M(T),Hcrit(T),T) where .Fis the 
free energy, /J, cr it (T) is the critical chemical potential and 
M is the solution of Eq. ([9]). As a result, at a fixed aver- 
age n two phases with uq(T) = n(M = 0, /i cr it(T), T) 
and fijj(T) = n(M(T), n crit (T),T) coexist as deter- 
mined by Eqations (|9ll0[) . The region of phase coexis- 
tence is shown in Fig.l. For comparison with the MF 
solutions we performed Monte-Carlo simulations of the 
model Eq.(3) in absence of the Coulomb forces. Due to 
strong fluctuations in 2D the critical temperature deter- 
mined from MC simulations is reduced by factor of ~ 2 
in comparison to the MF result. 



IV. COULOMB FRUSTRATED FIRST ORDER 
PHASE TRANSITION 

Let us now consider the role of the Coulomb repul- 
sion. The area under the T cr i t (n) in Fig.l is the area of 
phase coexistence. If we fix the temperature the two 
phases with the bulk concentrations hq and Um will 
have volume fractions 1 — x and x respectively, where 
x = (n — no)/(riM — no). Since the system is globally 
electroneutral, the phases with no and Um are charged. 
However, to break electroneutrality requires a large in- 
crease of the Coulomb energy. As a consequence growth 
of charged regions with two different charge densities is 
blocked by the Coulomb forces. 

In the literature there are a few examples of introduc- 
tion of c harg ing effects in the problem of phase separa- 
tion [H, 134 l35l |3H | . There are several different pos- 
sibilities to include long-range Coulomb forces in the 
model. Muratov[3(| proposed that the order parameter 
is a charged scalar and the charge density is proportional 
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to the order parameter. This situation is similar to the 
problem of a charged Bose gas in magnetic field consid- 
ered in [5§]. Similar situation is considered in Ref.(3f| 
where the free energy has two distinct minima as a func- 
tion of the density and gradient terms in the free energy 
are replaced by the surface tension. Jamei, Kivelson and 
Spivak [34| considered the case with a scalar order pa- 
rameter where the charge density is coupled to the order 
parameter as an external field. 

In our case, symmetry allows coupling of the charge 
density with square of the order parameter only. Let us 
consider the classical free energy density corresponding 
to the first order phase transition: 



Fx = ((t 1) + ( V 2 - 1) V- 



(11) 



Here t = (T — T c ) / (To — T c ) is the dimensionless tempera- 
ture. At t = 4/3 (T = T + (To - T c ) /3) a nontrivial min- 
imum in the free energy appears. At t = 1 (T = T ) the 
first order phase transition occurs. Below t = 1 the triv- 
ial solution r\ = corresponds to the metastable phase. 
At t = (T = T c ) the trivial solution becomes unsta- 
ble. In order to study the case of the Coulomb frustrated 
phase transition we have to add coupling of the order pa- 
rameter to the local charge density. In our case the order 
parameter describes the sublattice orbital magnetization 
and therefore only square of the order parameter can be 
coupled to the local charge density p: 



F, 



coupl 



-ar] 2 p. 



(12) 



The proposed free e nerg y functional is similar to that 
proposed in the Ref . [3J| . In our case the charge plays a 
role of the local temperature, while in Ref.[34| there is a 
linear coupling of the charge to the order parameter and 
the charge density plays a role of the external field. 

The total free energy density should also contain the 
gradient term and the electrostatic energy: 

Fgrad + Fa = C(Vr)) 2 + \k[ P (v) - p] 

dv'[p(v')-p]/\r-r\. (13) 

Here we write p explicitly to take into account global 
electroneutrality. The total free energy Eqs. (|llll2ll3|) 
should be minimized at fixed t and p. 

Next, we proceed to show that the Coulomb term leads 
to phase separation in 2D. Minimization of F with re- 
spect to the charge density p(r) leads to the following 
equation: 



aV 2 DV 2 = AttK[p(v) - p]dS(z), 



(14) 



here we write explicitly that electrostatic field is 3D but 
the charge density p(r) is confined in the 2D plane (z = 0) 
and depends only on 2D vector r. To preserve correct 
dimensionality we introduce the layer thickness d. We 
believe that this condition is favorable for creating charge 
segregation because electrostatic field is not screened in 



the third direction. Solving this equation by applying 
the Fourrier transform and substituting the solution back 
into the free energy density we obtain: 



F = F\ — olt] p - 
a 2 



iir 2 Kd 



dr 



C(V V ) 2 - 
Vfa(r) 2 )Vfa(r') 2 ) 



(15) 



As a result the free energy functional is similar to the 
case of first order phase transition with a shifted critical 
temperature due to presence of the term arfp and with 
additional nonlocal gradient term. 

To demonstrate that the uniform solution has a higher 
free energy then a nonhomogeneous solution we make the 
Fourier transformation of the gradient term: 



F gr ad OC Cfc 2 |77 k | 2 



« 2 fc|(?7 2 )k 
AirKd 



(16) 



where r/^ an d (f? 2 )k are Fourier components of the or- 
der parameter and square of the order parameter respec- 
tively. If we assume that the solution is uniform i.e. 
r/o 7^ and (r] 2 )o 7^ small nonuniform corrections to 
the solution reduce the free energy at small k, where the 
second term dominates. 

The situation is different in 3D. Direct solution of the 
equation for the charge density leads to the local gradient 
term of higher order — g^(Vr;(r) 2 ) 2 . This term can also 
lead to instability and higher order expansion in gradient 
terms become important. 



V. MONTE CARLO SIMULATIONS 

To substantiate above arguments we performed Monte- 
Carlo (MC) simulations of the system described by the 
Hamiltonian Eq. (3) with and without the presence of the 
long-range CR. The simulations were performed on a 
square lattice with dimensions L x L sites with 10 < 
L < 100 at different dimensionless temperatures t — 
k B Te a a/e 2 . The short range potential w/(i) = V/(i)eoa/e 
was taken to be nonzero only for |i| < 2 and was there- 
fore specified by two parameters: t>j(l, 0) and vi(l, 1) . To 
further minimize the number of free parameters only the 
nearest neighbors attractive interaction potential 17 (1, 0) 
was taken to be nonzero in most cases presented here. 

We first performed MC simulations of the model at a 
constant chemical potential in the absence of CR. Due to 
presence of first order phase transition, the particle den- 
sity probability distribution P t . M (n) has two peaks when 
the chemical potential is near the critical value p C rit(t)- 
At p C rit(t) the two peaks have equal height correspond- 
ing to the densities of the two coexisting phases, no and 
um- A standard Metropolis algorithm [50l| in combination 
with simulated annealing 51| and histogram reweighting 
technique [5^1 gave reliable results only at higher temper- 
atures near maximum t cr it(n). To improve reliability at 
lower temperatures we used a variant of multicanonical 
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aproachJH^I adapted to uniformly sample states over the 
full range of densities (58j. n, at a constant dimension- 
less temperature t S i m and chemical potential jU s j m . At 
each temperature the final histogram aquisition run in- 
volved at least 10 6 MC pseudospin flips per site. The 
density probability histograms Pt sim ,^{n) for several val- 
ues of the chemical potential n close to the simulation 
chemical potential a S i m were then calculated at each t S i m 
by reweighting [56] (see Fig. 2). From the histogram 
with equal peak heights the densities of the two coex- 
isting phases Uq and um were then determined at given 

In simulations at constant n one MC step consisted 
from a single update per each particle, where the trial 
move consisted from setting S z = at the site with 
nonzero Qi and S z — ±1 at a randomly selected site with 
zero Q{. A standard Metropolis algorithm [5C| in com- 
bination with simulated annealing [5l| was used in this 
case. A typical simulated annealing run consisted from 
a sequence of MC simulations at different temperatures. 
At each temperature equilibration phase consisting from 
10 3 — 10 6 MC steps was first executed followed by the av- 
eraging phase consisting from the same or greater num- 
ber of MC steps. Observables were measured after each 
averaging MC step during the averaging phase only. 

At constant n in absence of the CR global phase sepa- 
ration below t cr it (n) occurs in the form of a large cluster 
with M 0. To detect onset of clustering we measure the 
nearest neighbor density correlation function (CF) g p L = 
4»(i-»)L^ E|m|=i (Ei (Qi+m ~ n) (Qi - n)) L , where Q L 
represents the MC average. We define the temperature 
t c i(n) at which g p ^ rises to 50% of its low temperature 
value (see Fig. lb) as the characteristic crossover tem- 
perature related to the formation of clusters. 

In Fig. la we show the results of MC simulations 
in absence of the Coulomb repulsion. We find that for 
n > 0.4 the boundary conditions strongly affect the 
t C rit(n) line calculated at the constant chemical potential. 
When we use open boundary conditions (OBC) t crit (n) 
is strongly supressed above n 0.4 in comparison to the 
result obtained with the periodic boundary conditions 
(PBC). At constant n, on the other hand, the influence 
of the boundary conditions on t c i(n) is less pronounced. 
The t c i (n) calculated with both types of boundary condi- 
tions closely follow the t cr i t (n) line calculated with PBC. 
Above n ss 0.6 t c i(n) for OBC is only slightly higer than 
for PBC. We attribute insensitivity of t c i(n) to bound- 
ary conditions at fixed n to sensitivity of the correlation 
function to the short range correlations which are less 
sensitive to boundary conditions. 

Next we analyze the model in the presence of the long- 
range CR at constant n. In Fig. 3 we show a typical tem- 
perature dependence of the average energy per particle 
estimator (eMc)^ for different system sizes L in presence 
of the long range CR using OBC. Error bars represent the 

standard deviation o eMC — J (e\ lc ) L — (£mc)\- The 
average energy monotonously drops with decreasing tem- 
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FIG. 1: (a) The phase diagram of the model in absence of 
the Coulomb repulsion. The dashed line represents the mean 
field (MF) solution. The full sqaures (■) represent the t cr it(n) 
line calculated for the periodic boundary conditions by means 
of the multicanonical MC algorithm with the system of size 
L = 40. For comparison the t cr it{n) line is shown (A) for open 
boundary conditions. The t c i(n) lines for periodic (□) and 
open (0) boundary conditions are also shown, (b) The 
dependence of the nearest neighbor density correlation 
function, g p L , on temperature in absence of the Coulomb 
repulsion. The definition of t c i (n) is indicated by arrow. 
The numerical errorbars are of the order of symbol sizes 



perature. The drop is more pronounced in the temper- 
ature interval ~ 0.5 > t >~ 0.1 in which clusters start 
to form. Below t ~ 0.1 the clusters are partially or- 
dered. The temperature dependence of (ejvfc)i ^ s virtu- 
ally identical for all L (We should note that the curves 
are vertically shifted by 0.1 for clarity.) indicating that 
the boundary effects on (eMc) l are negligible even for 
the smalest system sizes. 

To check reliability of our simulations we analysed MC 
update dynamics by calculating the autocorrelation func- 
tion of energy fluctuations, 



aeh{ruc) 



\ 



K 



(eMc(i + t~mc) 



(e M c) L ){e MC [i) - (eMc) L ), (17) 



where euc (*) represents the energy per site at i-th MC 
step and tmc represents the MC time. A typical time 
dependence of gehijMc) is shown in the inset in Fig. 
4. The autocorrelation function drops with the char- 
acteristic MC relaxation time tr. 1/tr displays Arhe- 
nius temperature dependence (see Fig. 4) down to the 
temperature where clusters start to order. Below this 
temperature tr behaves more erraticaly. The activation 
energy strongly depends on the magnitude of the short 
range potential uj(l, 0). 

In the temperature region where clusters partially or- 
der the heat capacity cl = d (cmc) l displays the 
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FIG. 2: Histograms of the density probability distribution 
Pt,fi(n) at the chemical potential near fj, cr it(t) obtained by 
multicanonical MC simualtion in absence of Coulomb repul- 
sion The values of the coexisting densities no and um at the 
given temperature are indicated by arrows. Note the logarit- 
mic scale. 




FIG. 4: The characteristic MC relaxation time tr as a func- 
tion of temperature for different values of 0). Thin 
lines represent the Arrhenius fits. The inset shows the au- 
tocorrelation function <? e 30 (i~mc) at a few temperatures for 
n; (1,0) = — 1. For convenience tr is defined as a value of 
tmc where g e h{TMc) = 0.25. 
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FIG. 3: A typical temperature dependence of the average en- 
ergy per particle estimator (cmc) l for different system sizes 
L. Insets show snapshots of particle distribution at different 
temperatures, where darker and brighter shades of grey repre- 
sent Si — 1 and —1 respectively. Curves are vertically shifted 
for 0.1, error bars represent cr ejMC . 



peak at t co (n) (see Fig. 5b). The peak displays no scaling 
with L indicating that no long range ordering of clusters 
appears. Inspection of the particle distribution snapshots 
at low temperatures 33] reveals that finite size domains 
form (see Fig. 6). Within the domains the clusters are 
ordered. The domain wall dynamics seems to be much 
slower than our MC simulation timescale preventing do- 
mains to grow. The effective L is therefore limited by the 
domain size. This explains the absence of the scaling and 
clear evidence for a phase transition near t co (n). From 
the simulations is therefore not clear whether the absence 
of complete cluster ordering is due to the finitenes of the 



MC simulation or it is also due to the glassy form of 
the free energy landscape. The square shape of the sam- 
ple my frustrate the cluster orders with nontetragonal 
symmetries, while practically achievable number of MC 
steps per temperature step warrant reliable MC averages 
only above the temperature which is of the same order 
ast co (n). The cluster ordering temperature t co {n) which 
is the lowest energy scale at all densities is only weakly 
n-dependent between 0.1 < n < 0.9. 

We now focus on the short range potential shape which 
promotes formation of stripes[33|]. We set u;(l,0) = — 1 
and vi(l, 1) = and study the dependence of clustering 
on particle density. To detect clustering we again use the 
nearest neighbour CF. In Fig. 5c we plot a typical near- 
est neighbor CF, g p 4o(l, 0), as a function of temperature. 
At high temperatures t » |u;(l,0)| CF is slightly neg- 
ative due to the long range CR. When the temperature 
decreases CF becomes positive and further rises with the 
decreasing temperature. No saturation of CF as in the 
case of absence of the CR forces is observed with the de- 
creasing temperature (see Fig. lb). Again we define the 
temperature at which CF rises to 50% of its low temper- 
ature value as the characteristic crossover temperature, 
t c i(n), related to the formation of clusters. The depen- 
dence of t c i(n) on the particle density is shown in Fig. 
5a for different boundary conditions. While in absence 
of the long-range CR t c i(n) closley follows the t cr it(n) 
line (Fig la), suppression of clustering by the CR forces 
results in a significant decrease of t c i(n). 

Different boundary conditions influence t c i(n) only for 
densities above n > 0.5. In this region the particles that 
form clusters are holes (Qi — 0) in the background of 
pseudospins (Qi = 1). The open boundary conditions are 
effectively a perimeter formed form holes which attracts 
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FIG. 5: (a) The phase diagram of the model in presence of the 
long range Coulomb repulsion. Open circles (o) and full cir- 
cles (•) represent the t c i(n) line for periodic and open bound- 
ary conditions respectively while dotted circles (0) represent 
the td (n) line for periodic boundary conditions in absence of 
the long range CR. The onset of clustering, to(n) is shown 
by open squares (□) and the cluster ordering temperature 
tco(n) by open triangles (A). The pseudospin (orbital) or- 
dering temperature is shown by full stars (~k). Note the loga- 
rithmic scale. The error-bars are of the order of symbol sizes 
or smaller. 
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n 

FIG. 6: The density dependence of the cluster-size distribu- 
tion function Xj for a few smallest cluster sizes as a function of 
the average density at the temperature t = 0.14. The regions 
of densities where pairs prevail are shadowed. 



holes and by pinning enhances hole clustering resulting 
in an increase of t c i(n) for OBC. 

In addition for n > 0.5 and our choice of vi(i) the 
pseudospin background ferromagnetically orders at ts (n) 
which increases with increasing density as shown in Fig. 
5a. The pseudospin ordering temperature is significantly 
higher than t c i(n). Despite this the particle- hole sym- 
metry of the t c i(n) line is absent. The absence of the 
the particle-hole symmetry is a consequence of different 
entropy contributions of doubly degenerate particle level 



(Six == ±1 for Qi — 1) an d singly degenerate hole level 
(Si, = for Qi = 0). 

The t c i(n) line does not appear smooth. There are 
clear dips at n w 0.14, n = 0.5 and n w 0.86. With 
increasing density the ground state of the system appar- 
ently goes through a series of crossovers related to the 
most probable cluster sizes as shown in Fig. 6. While 
the dip at half filling clearly corresponds to commensu- 
rate ordering of stripes the other two dips approximately 
correspond to the densities at which clusters of size four 
start to replace pairs (see Fig. 6.). There is no obvious 
commensuration to underlying lattice at these densities. 
At densities at which larger clusters start to replace fours 
no comparable anomaly is observed in the t c i(n) line. 

Despite the presence of the CR forces some clusters 
already start to form at temperatures higher than t c i (n) . 
We can estimate the upper limit for the onset of cluster 
formation by the temperature, to{n), at which g p i,(l,0) 
crosses 0. It is interesting that to(n) almost coincides 
with the t C rit(n) line (see Fig. 5a) below n < 0.4 while at 
higher densities the onset of clustering appears at much 
higher temperatures. In the region 0.5 < n < 0.75 the 
onset of clustering is higher in temperature than the pseu- 
dospin ordering temperature ts(n) while above n « 0.75 
the pseudospin ordering represents the highest energy 
scale. 

To get further insight in the cluster formation we mea- 
sured the cluster-size distribution function. In Fig. 6 
we show the low temperature density dependency of 
the cluster-size distribution function, xj = N p (j)/(nL 2 ), 
where N p (j) is the number of particles for n < 0.5 or 
holes for n > 0.5 in clusters of size j. At the highest 
temperature Xj is close to the distribution expected for 
the random ordering. When the temperature decreases 
the number of larger clusters starts to increase at the 
expense of the single particle number. [33] Further down 
in temperature depending on the average density n clus- 
ters of a certain size start to prevail at the expense of all 
other sizes. Depending on the particle density prevailing 
clusters can be pairs up to n « 0.14, quadruples up to 
n w 0.3 etc.. The situation is qualitatively symmetrical 
for the clusters formed by holes at n > 0.5. We should 
note that for the given i?;(l,0) the system prefers clus- 
ters with an even number of particles, however different 
vi(i,j) might lead to the preference for an odd number 
of particles in a cluster. 

It should also be emphasized that the preference to cer- 
tain cluster sizes becomes clearly apparent only at tem- 
peratures lower then td(n), however the transition is not 
abrupt but gradual with decreasing temperature. This is 
seen also from gradual increase of the average cluster size 
with decreasing temperature shown in Fig. 7. Around 
half filling the average cluster size starts to diverge at 
low temperatures indicating formation of long stripe-like 
objects (see insets in Fig. 6) and proximity of the perco- 
lation. 
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FIG. 7: TThe temperature dependence of the average cluster 
size for different particle densities n below half filling (particle 
clusters) (a) and above half filling (hole clusters) (b). 



VI. CONCLUSIONS 

We presented the results of extensive investigation of 
the ordering of charged Jahn- Teller polarons as a function 
of doping and temperature. We consider charged parti- 
cles on a 2D square lattice subject to only the long-range 
Coulomb interaction and an anisotropic Jahn- Teller (JT) 
deformation. 

We prove that without the long-range Coulomb repul- 
sion the system is unstable with respect to the first or- 
der phase transition below the density dependent critical 
temperature. This was demonstrated by the solution of 
the mean field equation as well as by direct Monte Carlo 



simulations. It was shown that this result does not de- 
pend on the type of boundary conditions and the error 
due to finite size effect is estimated. 

In the presence of the Coulomb repulsion the global 
phase separation becomes unfavorable and the system 
shows a mesoscopic phase separation, where the size of 
the charged regions is determined by the competition 
between the ordering energy and the Coulomb energy. 
The phenomenological theory of this effect was formu- 
lated where the square of the order parameter is coupled 
with the charge density. The charge density plays the 
role of the local temperature. This type of coupling is 
more general in comparison with the models where the 
charge plays the role of an external held. 

Using Monte-Carlo (MC) simulations we showed that 
below a characteristic clustering temperature the system 
forms many different mesoscopic textures, such as clus- 
ters and stripes, depending only on the magnitude of the 
Coulomb repulsion compared to the anisotropic lattice 
attraction and the density of charged particles. Below 
the clustering temperature the system goes through a 
series of crossovers between phases with different meso- 
scopic textures when the particle density is increased. 
The low temperature part of the phase diagram is rather 
symmetric with respect to half filling. However, above 
half doping another high temperature scale appears cor- 
responding to orbital ordering of the particles. Surpris- 
ingly, a feature arising from the anisotropy introduced by 
the Jahn- Teller interaction is that objects with an even 
number of particles more stable than those with an odd 
number of particles. Such a behaviour could have signif- 
icant implications for superconductivity when tunnelling 
is included [13]. 



[1] E.Dagotto, Rep.Mod.Phys. 66, 763 (1994) 

[2] J. Jaklic and P. Prelovsek, Adv. Phys. 49, 1 (2000). 

[3] R.J.McQueeney et al, Phys. Rev. Lett, 82, 628 (1999). 

[4] A.Bianconi et al Phys.Rev.Lett 76 3412 (1996), Bozin et 
al., Phys.Rev.Lett. 84, 5856 (2000). 

[5] S.H.Pan et al., Nature 413, 282 (2001), McElroy et al., 
Nature 422, 592 (2003). 

[6] J.Demsar et al, Phys. Rev. Lett. 82, 4918 (1999), see 
also review by D.Mihailovic and V.V.Kabanov in "Su- 
perconductivity", ACS series on Structure and Bond- 
i ng. Eds. A .Bussmann-Holder and K.A.Miiller, (2004). 
, |cond-mat/040720"4| 

[7] see also: D.Mihailovic and K.A.Miiller, in High T c Su- 
perconductivity 1996: Ten years after the discovery. Eds. 
E.Kaldis, E.Liarokapis and K.A.Miiller, NATO ASI, Ser. 
E. Vol. 343 (Kluwer, 1997) p. 243. 

[8] J.M. De Teresa, M.R. Ibarra, P.A. Algarabel, C. Ritter, 
C. Marquina, J. Blasco, J. Garcia, A. del Moral, and Z. 
Arnold, Nature (London) 386, 256 (1997). 

[9] G. Allodi, R. De Renzi, G. Guidi, F. Licci, and M. W. 
Pieper, Phys. Rev. B 56, 6036 (1997). 
[10] M. Uehara, S. Mori, C. H. Chen, and S.-W. Cheong, 



Nature (London) 399, 560 (1999). 
[11] M. Fath, S. Freisem, AA. Menovsky.Y. Tomioka, J Aarts, 

JA Mydosh, Science 285, 1540 (1999). 
[12] G. Papavassiliou, M. Fardis, M. Belesi, T. G. Maris, G. 

Kallias, M. Pissas, D. Niarchos, C. Dimitropoulos, and 

J. Dolinsek, Phys. Rev. Lett. 84, 761 (2000). 
[13] Y. Kawasaki, T. Minami, Y. Kishimoto, T. Ohno, K. 

Zenmyo and H. Kubo T. Nakajima and Y. Ueda, Phys. 

Rev. Lett. 96, 037202 (2006). 
[14] P. L. Kuhns, M. J. R. Hoch, W. G. Moulton, A. P. Reyes, 

J. Wu, and C. Leighton, Phys. Rev. Lett. 91, 127202 

(2003). 

[15] M.H. Sage, G.R. Blake, G.J. Nieuwenhuys, T. T. M. Pal- 

stra, Phys. Rev. Lett. 96, 036401 (2006). 
[16] S-W. Cheong, H. Y. Hwang, C. H. Chen, B. Batlogg, L. 

W. Rupp, Jr., and S. A. Carter, Phys. Rev. B 49, 7088 - 

7091 (1994). 

[17] E. Dagotto, T. Hotta, A. Moreo, Physics Reports, 344, 
1 (2001). 

[18] J. Zaanen, O Gunnarsson, Phys. Rev. B, 40, 7391 (1989). 
[19] V.J. Emery, S. Kivelson and O.Zachar, Phys.Rev.B 56, 
6120 (1997). 



9 



[20] L.P. Gorkov, A.V. Sokol, Pisma ZhETF, 46, 333 (1987); 

L.P. Gorkov, J. Supercond., 14, 365, (2001). 
[21] J. Vitins, P. Wachter, Phys. Rev. B12, 3829, (1975). 
[22] Y. Shapira, S. Foner, N. Oliveira, T. Reed, Phys. Rev. 

B, 5, 2674 (1972). 
[23] Y. Shapira, S. Foncr, N. Oliveira, Phys. Rev. B, 10, 4765 

(1974). 

[24] V.J.Emery, S.A.Kivelson, H. Lin, Phys. Rev. Lett. 64, 
475 (1990) 

[25] V.J.Emery, S.A.Kivelson, Physica C, 209, 597 (1993). 
[26] G. Uhrig, R. Vlaming Phys. Rev. Lett., 71, 371 (1993). 
[27] A. Angelucci, S. Sorella Phys. Rev. B, 47, 8858 (1993). 
[28] A. Singh, Z. Tesanovic, J. Kim, Phys. Rev. B, 44, 7757 
(1991). 

[29] Y. Bang, G. Kotliar, C. Castellani, et al, Phys. Rev . B, 

43, 13724 (1991). 
[30] A. Moreo, D. Scalapino, E. Dagotto, Phys. Rev. B, 

textbf43, 11442 (1991) 
[31] U.Low, V.J.Emery, K.Fabricious, S.A.Kivelson, Phys. 

Rev.Lett. 72,1918 (1994). 
[32] J. Lorenzana, C. Castellani, and C. Di Castro, Phys. Rev. 

B, 64, 235127 (2001); Europhys. Lett. 57, 704 (2002). 
[33] T. Mertelj, V. V. Kabanov, D. Mihailovic, Phys. Rev. 

Lett. 94, 147003 (2005). 
[34] R. Jamei, S. Kivelson, and B. Spivak, Phys. Rev. Lett., 

94, 056805 (2005). 
[35] C Ortix, J. Lorenzana, C. Di Castro, Phys. Rev. B, 73, 

245117 (2006). 
[36] C.B. Muratov, Phys. Rev. E, 66, 066108 (2002). 
[37] D. Mihailovic, V.V.Kabanov and K.A.Muller, Europhys. 

Lett. 57, 254 (2002). 
[38] A.S. Alexandrov, A.M. Bratkovsky, and V.V. Kabanov, 

Phys. Rev. Lett., 96, 117003 (2006). 
[39] F.V. Kusmartsev, Phys. Rev. Lett., 84, 530, (2000), ibid 

84, 5026 (2000). 
[40] A.S. Alexandrov, V.V. Kabanov Pisma ZhETF, 72, 825 

(2000) (JETP letters, 72, 569 (2000)). 



[41] G.M. Zhao, K. Conder, H. Keller, K.A. Muller, Nature, 

381, 676 (1996). 
[42] G.M. Zhao, M.B. Hunt, H. Keller, K.A. Muller, Nature, 

385, 236 (1997). 
[43] D. Mihailovic, V.V. Kabanov, Phys. Rev., B, 63, 054505, 

(2001); V.V. Kabanov, D. Mihailovic, Phys. Rev., B, 

65, 212508, (2002), V.V.Kabanov and D. Mihailovic, 

J.Superc. 13, 959 (2000). 
[44] D.I. Khomskii, K.I. Kugel, Europhys. Lett. 55, 208 

(2001); Phys. Rev. B, 67, 134401 (2003). 
[45] M.B. Eremin, A.Yu. Zavidonov, B.I. Kochelaev, ZhETF, 

90, 537 (1986). 

[46] S.R. Shenoy, T. Lookman, A. Saxena, A.R. Bishop, Phys. 

Rev. B, 60, R12537 (1999); T. Lookman, et al, Phys. 

Rev. B, 67, 024114 (2003). ;K.H. Ahn, T. Lookman, A.R. 

Bishop, Nature 428, 401 (2004). 
[47] J. Lajzerovicz, J. Sivardiere, Phys. Rev. All, 2079 

(1975). 

[48] J. Sivardiere, J. Lajzerovicz, Phys. Rev. All, 2090 (1975) 
[49] V.V. Kabanov, A.S. Alexandrov, Phys. Rev. B, 71, 
132511 (2005). 

[50] N. Metropolis, et al. J. Chem. Phys. 21, 1087 (1953). 
[51] S. Kirkpatrick, CD. Gelatt and M.P. Vecchi, Science 220 
(1983) 671-680. 

[52] E.M.Lifshitz and L.P.Pitaevski, Physical Kinetics, ch.12 

(Buttcrworth-Heinemann, 1980). 
[53] T.Mertelj et al. (to be published). 
[54] D.Reagor et al., Phys.Rev.Lett. 62, 2048 (1989). 
[55] A.S.Alexandrov and N.F.Mott "Polarons and Bipo- 

larons", (World Scientific, 1995). 
[56] A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 

61, 2635 (1988). 
[57] B. A. Berg, T. Neuhaus, Physics Letters B 267, 249 

(1991). 

[58] G. Orkoulas and A. Z. Panagiotopoulos, J. Chem. Phys. 
110, 1581 (1999). 



